/*
Copyright (c) 2009-2010, Warwick Warp Limited
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of the Warwick Warp Limited nor the
      names of its contributors may be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL WARWICK WARP LIMITED BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#ifndef WARRAY_WARRAY_HELP_H_INCLUDED
#define WARRAY_WARRAY_HELP_H_INCLUDED

#include <complex>
#include <cmath>

#include "wassert.h"
#include "warray_nd.h"
#include "warray_sys.h"

namespace warp {
namespace detail {

template<typename LT,typename RT,warray_types::size_type N,typename Lptr,typename Rptr>
void bounds_check(const const_arraynd<LT,N,Lptr>& l, const const_arraynd<RT,N,Rptr>& r)
{
#ifdef BOUNDS_CHECK
	warp_assert_equal(l.dims(), r.dims());
#endif
}

template<typename T1,typename T2,typename T3,warray_types::size_type N,typename P1,typename P2,typename P3>
void bounds_check(const const_arraynd<T1,N,P1>& x, const const_arraynd<T2,N,P2>& y, const const_arraynd<T3,N,P3>& z)
{
#ifdef BOUNDS_CHECK
	warp_assert_equal(x.dims(), y.dims());
	warp_assert_equal(x.dims(), z.dims());
#endif
}

template<typename T>
void alias_check(const T*x, int incx, size_t nx, const T*y, int incy, size_t ny)
{
#ifdef BOUNDS_CHECK
	const T* xmin = std::min(x,x+incx*nx);
	const T* xmax = std::max(x,x+incx*nx);
	const T* ymin = std::min(y,y+incy*ny);
	const T* ymax = std::max(y,y+incy*ny);
	warp_assert_true(xmax <= ymin || ymax <= xmin);
#endif
}

template<typename T>
class multiplies_conj : public std::binary_function< std::complex<T>, std::complex<T>, std::complex<T> >
{
public:
	std::complex<T> operator()(const std::complex<T>& l, const std::complex<T>& r) const
	{
		return l*std::conj(r);
	}
};

template<typename T>
class inv : public std::unary_function<T,T>
{
public:
	T operator()(const T& v) const
	{
		return T(1)/v;
	}
};

template<typename T>
class invsqrt : public std::unary_function<T,T>
{
public:
	T operator()(const T& v) const
	{
		return T(1)/std::sqrt(v);
	}
};

template<typename T>
class hypot : public std::binary_function<T,T,T>
{
public:
	T operator()(const T& l, const T& r) const
	{
		if(l<r)
		{
			T c = r/l;
			return std::abs(l)*std::sqrt(T(1) + c*c);
		}
		else
		{
			T c = l/r;
			return std::abs(r)*std::sqrt(T(1) + c*c);
		}
	}
};

template<typename T>
class cis : public std::unary_function<T,std::complex<T> >
{
public:
	std::complex<T> operator()(const T& v) const
	{
		return std::complex<T>(std::cos(v), std::sin(v));
	}
};

// We need these because the overloads for complex types are templated
template<typename T>
struct choose_sin : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::sin(v); }
};

template<typename T>
struct choose_cos : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::cos(v); }
};

template<typename T>
struct choose_tan : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::tan(v); }
};

template<typename T>
struct choose_asin : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::asin(v); }
};

template<typename T>
struct choose_asin< std::complex<T> > : public std::unary_function<std::complex<T>,std::complex<T> >
{
	std::complex<T> operator()(const std::complex<T>& v) const { return boost::math::asin(v); }
};

template<typename T>
struct choose_acos : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::acos(v); }
};

template<typename T>
struct choose_acos< std::complex<T> > : public std::unary_function<std::complex<T>,std::complex<T> >
{
	std::complex<T> operator()(const std::complex<T>& v) const { return boost::math::acos(v); }
};


template<typename T>
struct choose_atan : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::atan(v); }
};

template<typename T>
struct choose_atan< std::complex<T> > : public std::unary_function<std::complex<T>,std::complex<T> >
{
	std::complex<T> operator()(const std::complex<T>& v) const { return boost::math::atan(v); }
};

template<typename T>
struct choose_sinh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::sinh(v); }
};

template<typename T>
struct choose_cosh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::cosh(v); }
};

template<typename T>
struct choose_tanh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::tanh(v); }
};

template<typename T>
struct choose_asinh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return boost::math::asinh(v); }
};

template<typename T>
struct choose_acosh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return boost::math::acosh(v); }
};

template<typename T>
struct choose_atanh : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return boost::math::atanh(v); }
};

template<typename T>
struct choose_conj : public std::unary_function<std::complex<T>,std::complex<T> >
{
	std::complex<T> operator()(const std::complex<T>& v) const { return std::conj(v); }
};

template<typename T>
struct choose_abs : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::abs(v); }
};

template<typename T>
struct choose_abs< std::complex<T> > : public std::unary_function<T,std::complex<T> >
{
	T operator()(const std::complex<T>& v) const { return std::abs(v); }
};

template<typename T>
struct choose_sqrt : public std::unary_function<T,T>
{
	T operator()(const T& v) const { return std::sqrt(v); }
};

template<typename T>
struct choose_pow : public std::binary_function<T,T,T>
{
	T operator()(const T& x, const T& y) const { return std::pow(x,y); }
};

template<typename T>
struct choose_exp : public std::unary_function<T,T>
{
	T operator()(const T& x) const { return std::exp(x); }
};

template<typename T>
struct choose_log : public std::unary_function<T,T>
{
	T operator()(const T& x) const { return std::log(x); }
};

template<typename T>
struct choose_log10 : public std::unary_function<T,T>
{
	T operator()(const T& x) const { return std::log10(x); }
};

template<int A, size_t K>
struct const_pow
{
	enum { value = const_pow<A,K-1>::value*A };
};

template<int A>
struct const_pow<A,0>
{
	enum { value = 1 };
};

template<typename T> struct ptr_detect {};
template<typename T> struct ptr_detect<T*> { typedef T type; };
template<typename T> struct ptr_detect<const T*> { typedef const T type; };

template<typename Tptr>
class ptr_with_stride
{
public:
	ptr_with_stride(Tptr x, int incx) : m_x(x), m_inc(incx) {}
	
	ptr_with_stride& operator--()
	{
		m_x -= m_inc;
		return *this;
	}
	
	ptr_with_stride operator--(int)
	{
		ptr_with_stride copy(*this);
		m_x -= m_inc;
		return copy;
	}
	
	ptr_with_stride& operator++()
	{
		m_x += m_inc;
		return *this;
	}
	
	ptr_with_stride operator++(int)
	{
		ptr_with_stride copy(*this);
		m_x += m_inc;
		return copy;
	}
	
	typename ptr_detect<Tptr>::type& operator*()
	{
		return *m_x;
	}
	
	typename ptr_detect<Tptr>::type& operator->()
	{
		return *m_x;
	}
	
	typename ptr_detect<Tptr>::type& operator[](ptrdiff_t index)
	{
		return m_x[index*m_inc];
	}
	
	ptr_with_stride& operator+=(ptrdiff_t offset)
	{
		m_x += offset*m_inc;
		return *this;
	}
	
	ptr_with_stride& operator-=(ptrdiff_t offset)
	{
		m_x -= offset*m_inc;
		return *this;
	}
	
	ptr_with_stride operator+(ptrdiff_t offset)
	{
		ptr_with_stride copy(*this);
		copy += offset;
		return copy;
	}
	
	ptr_with_stride operator-(ptrdiff_t offset)
	{
		ptr_with_stride copy(*this);
		copy -= offset;
		return copy;
	}
	
protected:
	Tptr m_x;
	int m_inc;
};

} // namespace
} // namespace

#endif
